Center for Turbulence Research 
Proceedings of the Summer Program 2004 


35 


Anisotropy of MHD turbulence at low magnetic 

Reynolds number 

By 0. Zikanov and A. Vorobev f, A. Thess \ P. A. Davidson and 

B. Knaepen || 


Turbulent fluctuations in MHD flows are known to become dimensionally anisotropic 
under the action of a sufficiently strong magnetic field. We consider the technologi- 
cally relevant case of low magnetic Reynolds number and apply the method of DNS of 
forced flow in a periodic box to generate velocity fields. The analysis based on differ- 
ent anisotropy characteristics shows that the dimensional anisotropy is virtually scale- 
independent. We also find that, except for the case of very strong magnetic field, the flow 
is componentally isotropic. Its kinetic energy is practically uniformly distributed among 
the velocity components. 


1. Introduction 

Magnetohydrodynamic (MHD) turbulent flows are ubiquitous in the universe, occur- 
ring in numerous astrophysical, geophysical, and technological applications. It is known 
that in the presence of a sufficiently strong magnetic field the turbulent fluctuations 
become anisotropic, which implies important consequences for the properties of the tur- 
bulence and possibly requires modification of numerical models. Specific manifestation 
of the anisotropy may vary but the principal mechanism is always the elongation of flow 
structures (turbulent eddies) along the lines of the magnetic field. 

The main difference between the kinds of MHD turbulence is due to different values 
of the magnetic Reynolds number 

Re™ = — , (1.1) 

V 

where r] = {aft o) _1 is the magnetic diffusivity, a and Ro being the electric conductivity of 
the liquid and magnetic permeability, and u , L are the typical velocity and length scales 
of the flow. If Re m > 1, there is a two-way coupling between fluctuations of magnetic field 
and velocity. This happens in astrophysical applications (stars, interstellar medium, etc.), 
where Re m 1, and in geophysical applications (Earth dynamo), where Re m is smaller 
but still significantly larger than 1. Discussions of anisotropy effects at large magnetic 
Reynolds number can be found, for example, in the recent review by Cho et al. (2002). 

The opposite case of Re m <C 1 occurs in a majority of technological processes, where 
a strong steady magnetic field is imposed on an electrically conducting liquid. Exam- 
ples of such applications include continuous casting of steel and aluminum, growth of 
semiconductor crystals, and lithium cooling blankets for fusion reactors. In this case, the 
low-Re m (so-called quasi-static) approximation can be applied (see, e.g., Moreau 1990 or 
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Davidson 2001). It can be assumed that the fluctuations of the magnetic field b associ- 
ated with fluid motions adjust instantaneously to the velocity fluctuations and that their 
amplitude is negligible in comparison with the amplitude of imposed magnetic field B. 
The rotational part of the Lorentz force reduces to the linear functional of the velocity 


F[v] = - 



<9 2 v 


(1.2) 


where p and u are the density and velocity of the fluid, A -1 is the reciprocal Laplace 
operator, and we assumed that the imposed magnetic field is uniform and purely vertical 
B = Be z . 

The flow transformation under the impact of force (1.2) has been actively studied in 
analytical (Moffatt 1967, Sommeria & Moreau 1982, Davidson 1997, 1999), experimen- 
tal (Votsish & Kolesnikov 1976, Alemany et al. 1979) and numerical (Schumann 1976, 
Zikanov & Thess 1998) works. The papers mentioned above represent only a fraction of 
the literature on the subject, more references being available therein. Far from the walls, 
the action of the magnetic field was identified as two-fold. First, the induced electric cur- 
rents result in additional dissipation of kinetic energy, the Joule (magnetic) dissipation. 
Second, the flow becomes anisotropic, its structures being elongated along the magnetic 
field lines. 

The reason for the anisotropy becomes especially transparent if one assumes that 
the flow is unbounded and uniform and uses the Fourier representation.! The Fourier 
transform of (1.2) is 


F[v] = -~ ( B ^ v(k, t) = — ^^-v(k, t) cos 2 9, (1.3) 

p k z p 

where k is the wavenumber vector and 9 is the angle between k and B. The rate of Joule 
dissipation of a Fourier mode with the wavenumber vector k is 


yu(k) = ^-(v(k, t) • v*(k, t)) cos 2 9, (1.4) 

P 

so the dissipation is anisotropic. It attains maximum for the Fourier modes with B || k 
and zero for modes with B ik, i.e., for modes independent of the ^-coordinate. The 
dissipation tends to eliminate the velocity gradients in the direction of B and elongate the 
flow structures in this direction. The limiting case is the two-dimensional state completely 
independent of the ^-coordinate. The Joule dissipation is equal to zero in this state. 

This picture of the flow transformation was first given by Moffatt (1967) and later 
beautifully illustrated by Sommeria & Moreau (1982), who introduced the image of ’Joule 
cone’ in the wavenumber space around the B-direction, in which the magnetic dissipation 
takes place. A remarkable feature of the picture is that, according to (1.4), the relative 
rate of the dissipation /r(k)/u 2 depends on the angle 9 but not on the wavenumber k. 
One can, thus, assume that the anisotropy would develop equally on all length scales of 
the flow. 

The situation, however, looks much more complicated if one takes into account the 
non-linearity of the Navier-Stokes equations and the resulting energy transfer between 
the modes and tendency to restoration of isotropy. The ratio between the Lorentz force 


f The process of development of anisotropy can be viewed from a different angle. In partic- 
ular, Davidson (1997) proposed a scenario based on the fact that the Lorentz force conserves 
momentum component parallel to the magnetic field. 
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and non-linear term is evaluated by the magnetic interaction parameter (Stuart number) 


N = 


oB 2 L 

pu 


(1.5) 


It is clear that the linearized picture of the flow development discussed above is true only 
in the limit of N 1. At finite N, one can expect a more complex scenario, probably 
with scale-dependent anisotropy. In particular, the analogy with stratified, rotating, or 
strained turbulence (see, e.g., Smith & Mourn 2000) suggests that smaller scales are 
more isotropic than large scales. Another way to arrive to the suggestion is to notice 
that, while the Joule damping time r = p/crB 2 is scale-independent, the eddy turnover 
time T = L/u is not. The magnitude of the scale-related Stuart number N = T/t can 
vary with the scale and so is the flow’s anisotropy. 

The question of anisotropy at different length scales is particularly important in the 
view of recent attempts to apply the traditional LES models to the low-Re m MHD 
turbulence (Knaepen & Moin 2004). The a-posteriori evaluation of the dynamic model 
(Germano et al. 1991) showed good accuracy for decaying turbulence at moderate hy- 
drodynamic Reynolds number and N < 10. This result is counter-intuitive since one 
expects LES models developed in assumption of local isotropy to perform poorly in the 
case of strongly anisotropic flow. One can not guarantee that equally good results will 
be obtained at higher Reynolds numbers. 

In this paper, we investigate the anisotropy of low-Re TO MHD turbulence using the DNS 
of forced flow in a box with periodic boundary conditions. The model and the numerical 
experiments are described in section 2. Various measures of anisotropy are discussed and 
evaluated in section 3. Possible implications for LES subgrid-scale modeling are discussed 
in section 4. Finally, summary and concluding remarks are provided in section 5. 


2. Model and numerical experiments 

We solve MHD equations for viscous, incompressible and electrically-conducting fluid 
in the quasi-static approximation. The Lorentz force term is given by (1.2). After apply- 
ing (Vx)x operation to eliminate pressure and taking Fourier transform, the governing 
equation become 


|jr(M) = -*v[ kx ( kx< i)]- 


uk 2 + 


oBl 

P 



(2.1) 


where v is the kinematic viscosity and q is the Fourier transform of the nonlinear term. 
The incompressibility condition used in derivation of (2.1) can be applied to recover the 
modified pressure field. A statistically homogeneous flow is calculated within a rectangu- 
lar box with periodical boundary conditions. Since we expect axial anisotropy of turbulent 
flow, and elongation of turbulent vortexes along the z-axis, the elongated box of dimen- 
sions 2n x 2n x 4n is used. The equation (2.1) is solved by the standard pseudo-spectral 
technique. 

In order to generate a statistically steady flow over long period of time, an artificial 
forcing is applied at large length scales. A constant amount of energy is added at each 
time-step. When equilibrium is reached the amount of energy added is on average equal 
to the energy dissipated. The external force 

H(k) = a(k)v(k) 


(2.2) 
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Figure 1 . (a), Evolution of total energy, (b), Evolution of rates of total viscous ( ) and 

magnetic ( ) dissipation. 


is applied to modes with 1.5 < k < 3.1. The time-dependent coefficients a are determined 
so that the net work by forcing is equal to the prescribed total (viscous and magnetic) 
rate of dissipation e 0 and the work is equally divided among the forced modes 

a(k) = — , (2.3) 

Nf orced (v(k) • v (k)) 

where Nf orce d is the number of forced modes and * stands for complex conjugate. 

Simulations are performed with numerical resolution 256 x 256 x 512. The param- 
eters eo = 0.5 and v = 2.2 ■ 10 -3 are chosen so as to guarantee the accuracy criterion 
kmaxV > 1-5, where k max is the maximum resolved wavenumber and r] is the Kolmogorov 
dissipation scale. The microscale Reynolds number is ReA ~ 94 in the non-magnetic run. 

The numerical experiments are staged as illustrated in figure 1. First, a developed 
turbulent flow is calculated starting with random, isotropic, and homogeneous field and 
continuing simulations without magnetic field for sufficiently long period. The complete- 
ness of the transitional period is judged by stabilization of the values of the total kinetic 
energy, viscous and magnetic dissipation rates defined as 

1 _d2 u2 

k k ' k 

(2.4) 

The flow field computed at the moment to = 9.72 is used as an initial condition for 
three simulations with different strength of the magnetic field. At t = to, the integral 
length scale is L = 0.73 and the rms turbulent velocity is u = 0.71, which gives the 
turbulent eddy turnover time of about 1. The strength of magnetic field in each run can 
be identified by the values of magnetic interaction parameter (1.5) at t = to- The cases 
with no magnetic field ( N 0 = 0), moderate magnetic field ( N 0 = 1), and strong magnetic 
field (No = 5) are consideredf. 

A comment must be made regarding the limitations of the traditional DNS approach 
that can affect the accuracy of our results. First, at moderate ReA, the artificial forcing 
can have some impact on a significant portion of the resolved length scales. Second 

f The flow evolution caused by magnetic field leads to significant increase of N during the 
runs (to about 1.2 in the second case and 10 in the third case). 
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Figure 2. Anisotropy coefficients G 1 ( ) and G2 ( ) given by (2.6). 


(perhaps less important for us here), the integral length scale is relatively large (about 
1/5 of the box size at t = t$). The periodic boundary conditions can, therefore, affect 
the flow behavior at largest scales. 

The evolution of total energy, viscous and magnetic dissipation is shown in figure 1. 
It can be seen that the period of adjustment after the introduction of magnetic field 
lasts several turnover times, after which the flow behavior is statistically steady. In the 
flows with magnetic field, the Joule dissipation is responsible for major part of the total 
dissipation rate eo = e(t) + p(t). 

Zikanov & Thess (2004) proposed using anisotropic Taylor microscales 


A< = 



d/2 


0- + Sij){{dv j /dx i ) 2 ) 


(2.5) 


as global measures of anisotropy of velocity gradients. Here, any velocity component v-i 
can be used and (• • •) stands for volume averaging. For isotropic turbulence, A^ are all 
statistically equal to the usual Taylor microscale A = [15 vu 2 /e] 1 / 2 . The ratio (A_l/A||) 2 of 
length scales measured in the directions transverse and parallel to the magnetic field is 
equal to often used anisotropy coefficients (see, e.g., Schumann 1976 or Zikanov & Thess 
1998) shown in figure 2 

({ dv 2 /dz) 2 ^ ({ dvzfdzf ^ 

7 ry j 2-^ » • (2-6) 

2 ({dv 2 ldyf) {(dv 3 /dyf) 


Figure 3 shows skewness and flatness of longitudinal derivatives of the velocity com- 
ponents 

(( dvi/dxi ) 3 \ (( dvi/dxi ) 4 \ 

S * = ~T d/i* * = 7 77- ( 2 - 7 ) 

{[dvi/dxi) J {(dvi/dxi) ) 

In the isotropic flow at N 0 = 0, all three components of skewness and flatness are approx- 
imately equal to 5 m 0.5 and F « 5.4. In the presence of magnetic field, the components 
measured across and along the magnetic field differ considerably. For skewness, the trans- 
verse components remain close to the isotropic value, although some drop can be seen in 
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a) t>) 




Figure 3. Skewness (a) and flatness (b) coefficients (2.7) measured across ((Si + S 2 )/2 and 

(Fi + F 2 )/2 shown as ) and along (53 and F 3 shown as ) magnetic field. For No — 0, 

solid lines show (Si + 52 + 53 )/ 3 and (Fi + F 2 + Fs)/3 . 



Figure 4. Modified pressure fields in developed flows. 


the case of N 0 = 5. The parallel component decreases significantly at N 0 = 1 and even 
more so at N 0 = 5. This clearly indicates suppression of nonlinear energy transfer and 
dissipation of strong longitudinal velocity gradients by the magnetic field. For flatness, 
we also observe significant reduction in the case of N 0 = 5, more pronounced for the 
transverse component. 

The internal structure of flow field is illustrated in figure 4, where we show snapshots of 
modified pressure field made in developed flows at different magnetic field. The tendency 
to anisotropy is clearly seen, although even at iVo = 5 the flow is far from approaching 
two-dimensional form. 

The last issue to be discussed in this section is that of componental anisotropy, which 
is understood as anisotropy of the Reynolds stress tensor or velocity field (Kassinos et al. 
2001). As the Joule dissipation directly affects only velocity gradients along the magnetic 
field lines, the componental anisotropy is a secondary effect, which existence and strength 
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a) Velocity field 


b) Vorticity field 



Figure 5. Invariant maps for componental anisotropy of velocity (a) and vorticity (b) fields. 
Squares, triangles, and circles are for IVo — 0, 1, and 5, correspondingly. 


is far from being obvious. The componental anisotropy of field v is traditionally evaluated 
with the help of the traceless tensor (Lumley & Newman 1977) 


_ { VjVj ) 

13 M 



(2.8) 


which has two nontrivial invariants 6 rf = bijbji and 6£ 3 = bijbjkbki- The magnitude of 
q represents the degree of anisotropy, with q = 0 for an isotropic flow, and the upper 
boundary q = (l/ 22 + 2£ 3 ) corresponding to a purely two-component case. The type 
of the anisotropy is determined by the sign of £. If £ > 0, the flow is dominated by 
the velocity component parallel to the axis of symmetry (so-called prolate axisymmetric 
flow). Negative values of £ mean domination of two perpendicular components (oblate 
axisymmetric flow). 

We calculated the invariants at several time moments separated by few eddy turnover 
times. The results presented on the invariant map in figure 5a show that the velocity fields 
at N 0 = 0 and N 0 = 1 are fairly isotropic. The flow with strong magnetic field at N 0 = 5 
demonstrates a degree of anisotropy. It is interesting that the velocity field migrates 
from being dominated by two velocity components perpendicular to the magnetic field 
to being dominated by one parallel component. This behavior can be associated with 
slow transformation of the large-scale structures of the flow. 

We also calculated tensor (2.8) and invariants q and £ for the vorticity field. The results 
shown in figure 5b are easy to explain. As the vertical gradients of velocity are eliminated 
in flows with higher N, the vorticity field becomes dominated by the vertical component. 


3. Anisotropy at different length scales 

In this section we analyze the calculated flow fields to determine how and whether 
at all the dimensional anisotropy varies with the length scale. We start with the power 
spectra shown in figure 6. The steepening of the energy spectrum and accompanying 
decrease of viscous dissipation are well known phenomena associated with the anisotropy 
(see, e.g., experiments by Alemany et al. 1979 or simulations by Zikanov & Thess 1998). 
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a) b) 





Figure 6. Spectra of (a) kinetic energy, (b) rate of viscous and (c) magnetic dissipation, and 
(d) spectra of anisotropy coefficient (3.1). 


More interesting is the fact that the slope of Joule dissipation curves in figure 6c closely 
follows the slope of kinetic energy curves. This behavior detected earlier by Zikanov & 
Thess (1998) and confirmed here for larger Reynolds number gives, at least, a partial 
answer to our question. The ratio 

_ 3 Tj n(k) _ 3E|k-u* _ 3D 33 (k) 

~ 2 E(k) £uu* 2 E(k) 

can be considered as a measure of dimensional anisotropy at the wavelength k. In (3.1), 
the sum is over all wavenumber vectors in the shell k— 1/2 < |k| < A; + 1/2, Tj = p/aB 2 
is the Joule dissipation time, and D 33 is the component of the dimensionality tensor 
considered by Kassinos et al. (2001). The scaling factor is chosen so as to provide g = 1 
in an isotropic flow. It can be seen in figure 6d that, outside of the forced region, g(k) 
varies only slightly with k both at moderate and strong magnetic field. Even this slight 
decrease can be attributed to the effect of forcing. The magnitude of g(k) is close to that 
of the global anisotropy coefficients G\ and G 3 shown in figure 2. 

The use of Fourier wavenumbers is arguably not the best way to consider the length 
scale decomposition of the flow properties. In the following we concentrate on analyzing 
the flow anisotropy in the physical space. First, figures 7 and 8 show the two-point 
velocity correlations along the direction perpendicular and parallel to the magnetic field. 
For the correlations in figure 7 we assume axial symmetry and calculate, for example, 
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Figure 7. Correlations of velocity field in the direction perpendicular to the magnetic field. 




Figure 8. Correlations of velocity field in the direction parallel to the magnetic field. 


Qxx = (u(x)u(x + re±)) / (u 2 ), where {• • •) now stands for space averaging in x and 
averaging over all horizontal directions of ej_. 

One can see in figures 7 and 8 that imposing a magnetic field leads to growth of corre- 
lations not only in parallel but also in perpendicular direction. This is in agreement with 
the development of larger coherent structures in this case (see figures 4). Interestingly, 
as can be seen from comparison of plots (a) and (b) in both figures, the vertical velocity 
component experience much smaller increase than the horizontal components. 

Figure 9 can be considered a physical-space analog of figure 6d. We calculate the 
second-order structure functions T ,;(r) = {(^(x+rej)— u^(x)) 2 ) and evaluate the anisotropy 
at the physical scale r as the ratio 

2T 3 

Ti + X 2 ' 

One can see that, again, the anisotropy is virtually scale-independent. 


(3.2) 
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Figure 9. Anisotropy coefficient R(r) given by (3.2). 


4. Anisotropy of stress and strain tensors 

One of the goals of this work is better understanding of properties of MHD flows 
related to the subgrid-scale turbulence modeling in LES. One can assume that the tra- 
ditional SGS closures developed under the assumption of local isotropy are not suitable 
for strongly anisotropic MHD flows at high N. The opposite conclusion was made by 
Knaepen & Moin (2004) who performed LES of decaying turbulence at moderate Re and 
N up to 10. They found that the dynamic Smagorinsky model (Germano et al. 1991) 
is not less accurate for the MHD flows than for isotropic flows without magnetic field. 
Whether such a good performance is specific to the dynamic model and whether it is 
related to the scale-independence of anisotropy discussed above requires further investi- 
gation. So far, we have analyzed the computed data by calculating root-mean-squares of 
components of the Reynolds stress and rate of strain tensors of filtered velocity field 


Sij 

T ij 2 T~kkj 


where = ViVj — ViVj 


and 


= 2 


dvj 

dxj 


+ 


dVj 

dxi 


(4.1) 


(4.2) 


Here, the — stands for Fourier cut-off filtering at a wavenumber k. The results are pre- 
sented as functions of k in figure 10. We assume the axial symmetry and take the averages 
of statistically equal quantities. It can be seen that the rate of strain tensor becomes in- 
creasingly anysotropic with growing N. There is a significant and consistent difference 
between the amplitudes of, say, 5 12 and S 13 . On the other hand, the components of the 
Reynolds stress tensor are virtually indistinguishable for any N and for any filter cut-off 
k. 


5. Conclusions 

We performed simulations of forced homogeneous turbulence in a low -Re m MHD flow. 
Three numerical experiments with equal constant energy input but different strength of 
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a) N o =0 b) N„=1 c) N 0 =5 



Figure 10. Spectra of root-mean squares of components of stress (top) and rate-of-strain (bot- 
tom) tensors of filtered velocity field as function of the filter cut-off k. Here n = ir'ii + r| 2 ) /2, 
72 = ( T 33 )) r 3 = (ti2 + T%1 ) /2, T 4 = ^Ti 2 3 + r|i + 723 -f T 32 ) /4. Similar notations used for 
rate-of-strain tensor. For IVo = 0, n and T 2 curves coalesce, as well as r 3 and T4. 

applied magnetic field were performed. The dimensional and componental anisotropy of 
the flow was analyzed both in the global sense and as a function of length scale. The 
main conclusions are as follows. 

• While the flow develops strong dimensional anisotropy under the action of the Joule 
dissipation, its componental anisotropy remains insignificant. The kinetic energy remains 
approximately equally distributed over the velocity components. 

• The analysis performed both in spectral and physical space showed that the dimen- 
sional anisotropy is virtually scale-independent. We can not exclude the possibility that 
different behavior can be observed at small scales of flows with much higher Reynolds 
numbers. On the other hand, the extent of the range of constant anisotropy was found 
to be quite significant in our simulations (about 2 orders of magnitude). We can consider 
the constancy a robust feature of low-i?e TO MHD turbulence. 
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